!******************************************************************************************
SUBROUTINE counter2_u_fe(gn_initial2,iax,a_initial2,b_initial2,i_default_global2,bpres,hres,gnres,wres)

USE GLOBAL

IMPLICIT NONE

REAL(8),INTENT(IN)::a_initial2,b_initial2,gn_initial2
INTEGER,INTENT(IN)::i_default_global2,iax
REAL(8),INTENT(OUT)::bpres,hres,gnres,wres
REAL(8)::pres,cnres,ctres,cres,plow,phgh,ctmax,gnval
REAL(8)::taures,ares,wwres,utilc,utilg ,a0,b0, aux0
REAL(8)::bpval,wval,tauval,slack,ucost
REAL(8)::p0,h0,q,b,udef
REAL(8)::ctval,cnval,cval,q_fun,penalty,fval,exp_v_next,taxbar
REAL(8)::q02(nb2),b2(nb2),qval,qbval,qbmin,qbmax,aux2(1),cont_value,old_valuemax
INTEGER(4)::i,jj,im,ibp,iaux,feas,error,niter2,iqbmin,bout,ftype,root
LOGICAL::flag00
INTEGER(4)::i_num_star,i_num_star2,index_optmax
real(8)::sign_new,sign_old,sign_two,ctemp2,difc
REAL(8)::bpresold,hresold,gnresold,wresold

INTEGER::i_num,index_opt
INTEGER,PARAMETER::num=300
REAL(8)::ptemp(num),bound1,bound2,new_value,old_value

EXTERNAL q_fun,cont_value

!OBS: If govt b.c. nonbinding, it's always optimal to reduce debt
!       ==> govt b.c. always binds

i_num_star2 = 0
root=0

b0 = b_initial2
a0 = DEXP(a_initial2)
gnval  = gn_initial2
cnval  = 1d0-gnval  

qbmin  = MINVAL(q1(iax,:)*(bgrid2-(1d0-lambda)*b0))
qbmax  = 0d0

DO ibp=1,nb2         !flip upside down (to find highest MINLOC), i.e. highest zero
    q02(ibp) = q1(iax,nb2-ibp+1)
    b2(ibp)  = bgrid2(nb2-ibp+1) -(1d0-lambda)*b0
ENDDO

aux2   = MINLOC(q02*b2)     !trough of laffer curve   
iqbmin = aux2(1)
iqbmin = nb2-iqbmin+1

ctmax = a0 -qbmin + payoff*b0 

!b.c. holds
ftype = 200
plow = wbar/alpha
phgh =(1d0-omegac)/omegac*(ctmax/cnval)**(1d0+mu)

old_value = AFunction(ftype,a0,gnval,a0,plow)
old_valuemax = 10d+33 

i_num_star=0
index_opt = 0
index_optmax = 0
difc = (phgh-plow)/(REAL(num)-1d0)

1212 DO i_num=1+i_num_star,num
    ftype = 200
    ptemp(i_num) = plow+(phgh-plow)*(REAL(i_num)-1d0)/(REAL(num)-1d0)     
    new_value = AFunction(ftype,a0,gnval,a0,ptemp(i_num))
    
    sign_new = SIGN(1d0,new_value)
    sign_old = SIGN(1d0,old_value)
    sign_two = SIGN(1d0,sign_new*sign_old)

    if (DABS(new_value) <= DABS(old_valuemax)) then
        index_optmax = i_num
        old_valuemax = new_value        
    end if      
    i_num_star = i_num
    if (sign_two.EQ.-1d0) then
        index_opt = i_num
        old_value = new_value
        GOTO 1515
    ENDIF
    old_value = new_value

ENDDO       
   
1515 IF (index_opt.EQ.0) THEN
    bound1=ptemp(MAX(index_optmax-1,1))
    bound2=ptemp(MIN(index_optmax+1,num))
ELSE
    bound1=ptemp(MAX(index_opt-1,1))
    bound2=ptemp(MIN(index_opt,num))   
ENDIF
    
p0=BrentRoots(ftype,a0,gnval,a0,bound1,bound2,1d-8,1000,fval,niter2,error)  

IF (error.EQ.1) THEN 
    p0 = 0d0
    IF (i_num_star.LT.num)  go to 1212
ELSEIF (error.EQ.2) THEN
    PRINT*,'Problem: Counter 2 FE: More Iterations!'
    PAUSE
    STOP
ENDIF
    
wval   = alpha*p0  
aux0   = p0*gnval-taxrate*(a0+in*p0)
tauval = p0*gnval-aux0
ctval  = a0+aux0

qbval  = payoff*b0-aux0 

IF ((qbval.LT.qbmin)) THEN
    bout = 1            !govt b.c. not satisfied
ELSE
    iaux   = nb2

    DO WHILE ((iaux.GT.iqbmin).AND.(q1(iax,iaux)*(bgrid2(iaux)-(1d0-lambda)*b0).GT.qbval))        !find highest index with min value
        iaux = iaux-1
    ENDDO

    bout = 0
    bpval  = LINEAR_INT(q1(iax,iaux:iaux+1)*(bgrid2(iaux:iaux+1)-(1d0-lambda)*b0),bgrid2(iaux:iaux+1),qbval) 

    qval   = qbval/(bpval-(1d0-lambda)*b0)
    
    IF (qval.LT.0d0) THEN
        PRINT*,'Problem price q = 0 in counter2 FE',gn_initial2
        bpres = 0d0
        hres = 0d0
        gnres= 0d0
        wres = huge_n
        RETURN
        PAUSE
        STOP
    ENDIF
   
ENDIF

taxbar = taxrate*(a0+in*p0)

flag00 = ((p0.LE.0d0).OR.(wval.LT.wbar).OR.(bout.EQ.1).OR.(ctval.LE.0d0).OR.(tauval.GT.taxbar)) 

IF (.NOT. flag00) THEN
    pres  = p0
    hres  = 1d0
    bpres = bpval
    ctres = ctval
    taures= tauval
    wwres = wval
    gnres = gnval
    feas  = 1d0
ELSE
    ctres   = 0d0    
    feas    = 0d0
    taures  = tauval
    hres    = 0d0
    pres    = -1d0
    wwres   = wval
    bpres   = 0d0
ENDIF

cnres = MAX(1d-8,hres**alpha-gnval)

IF (mu .NE. 0d0) THEN
    cres  = (omegac*ctres**(-mu)+ (1d0-omegac)*cnres**(-mu))**(-1d0/mu)
ELSE
    cres  = (ctres**(omegac))*(cnres**(1d0-omegac))
ENDIF

cres = cres - (cnres .LE. 1d-6)*(-1d6)
cres = cres - (ctres .LE. 1d-6)*(-1d6)

IF (sigma.EQ.1d0) THEN
    utilc = DLOG(MAX(cres,1d-8))
ELSE
    utilc = MAX(cres,1d-8)**(1d0-sigma)/(1d0-sigma) 
ENDIF

IF (sigmag.EQ.1d0) THEN
    utilg = DLOG(MAX(gnres,1d-8)) 
ELSE
    utilg = MAX(gnres,1d-8)**(1d0-sigmag)/(1d0-sigmag) 
ENDIF

!negative penalty
penalty = 1d3*((pres.LE.0d0)*(0d0-pres)+(wwres.LT.wbar)*(wbar-wwres))
penalty = penalty+1d9*((bpres.LT.bmin)*(bmin-bpres))+1d2*(bpres.LT.bmin)
penalty = penalty +1d3*((gnres.GT.1d0)*(gnres-1d0)+(gnres.LT.0d0)*(0d0-gnres))

ucost = defaultgrid(i_default_global2)*udef(a_initial2)

wres = (1d0-chi)*utilc+chi*utilg  + ucost   !deactivate default cost
wres = feas*wres + (1d0-feas)*huge_n +penalty       

exp_v_next = beta*cont_value(bpres,a_initial2,i_default_global2)
wres = wres + exp_v_next

!x_global1 = 1d0
!gn_global1 = gnres
!b_next_global1 = bpres
!counter2_u_fe = wres

IF (root.EQ.0) THEN
    bpresold = bpres
    hresold = hres
    gnresold = gnres
    wresold = wres
    
    root=1
ELSEIF (root.GT.0) THEN
    IF (wres.GT.wresold) THEN
        bpresold = bpres
        hresold = hres
        gnresold = gnres
        wresold = wres
    ENDIF
root=root+1

ENDIF
    
IF (i_num_star.LT.num)     GO TO 1212

IF (wresold.GT.wres) THEN
    bpres = bpresold
    hres = hresold
    gnres = gnresold
    wres = wresold
ENDIF
    
RETURN

end SUBROUTINE counter2_u_fe
!***************************************************************************
SUBROUTINE counter2_u_unemp(gn_initial2,iax,a_initial2,b_initial2,i_default_global2,bpres,hres,gnres,wres)

USE GLOBAL

IMPLICIT NONE

REAL(8),INTENT(IN)::a_initial2,b_initial2,gn_initial2
INTEGER,INTENT(IN)::i_default_global2,iax
REAL(8),INTENT(OUT)::bpres,hres,gnres,wres
INTEGER,PARAMETER::num=100
REAL(8)::pres,cnres,ctres,cres,a0,b0
REAL(8)::taures,ares,wwres,utilc,utilg,uval
REAL(8)::pp1,bpval,wval,tauval,aux0
REAL(8)::p0x,bpvalx,wvalx,tauvalx,uvalx,hvalx,hf1,ctvalx,cnvalx
REAL(8)::p0,h0,q,b,x_global2old,objective_u_unempold,lowh,hghh,hzero
REAL(8)::ctval,cnval,gnval,hval,cval,penalty,fval,ctmax
INTEGER(4)::i,j,jj,feas ,error,niter2,try,numval,ftype
LOGICAL::flag00,flag00x
REAL(8)::q02(nb2),b2(nb2),qval,qbval,qbmin,qbmax,aux2(1),fvalue
REAL(8)::plow,phgh,taxbar,q_fun,aahat,aa1    
REAL(8)::qvalx,qres,ucost,udef,exp_v_next,cont_value,ptemp(num),old_value,new_value,ftol,xval(2),xi(2,2)  
real(8)::sign_new,sign_old,sign_two,ctemp2,difc,old_valuemax,bound1,bound2
INTEGER(4)::i_num,index_opt,i_num_star,i_num_star2,index_optmax,root
INTEGER(4)::im,ibp,iaux,iqbmin,bout,ibx,ITERF
EXTERNAL cont_value,q_fun

i_num_star2 = 0
root=0

b0 = b_initial2
a0 = DEXP(a_initial2)
gnval  = gn_initial2

qbmin  = MINVAL(q1(iax,:)*(bgrid2-(1d0-lambda)*b0))
qbmax  = 0d0

try = 1
lowh = (gnval+ 1d-4 )**(1d0/alpha)   !only one positive root: hf0 s.t. f'(hf0)=0 is always negative
!hf1  = (x_glob/kw/(wbar/alpha*gnval**(1d0-alpha))+gnval*(1d0+1d0/kw))**(1d0/alpha)   
hf1  = 1d0         
hghh = MIN(1d0,hf1) 

ctmax = a0-qbmin + payoff*b0
plow = wbar/alpha*(gnval**((1d0-alpha)/alpha))
phgh = wbar/alpha

ftype=35
old_value = 10d+33 
old_value = AFunction(ftype,a0,gnval,gnval,plow)
       
i_num_star=0
index_opt = 0
index_optmax = 0
difc = (phgh-plow)/(REAL(num)-1d0)

1313 DO i_num=1+i_num_star,num
    ftype=35
    ptemp(i_num) = plow+(phgh-plow)*(REAL(i_num)-1d0)/(REAL(num)-1d0)     
    new_value = AFunction(ftype,a0,gnval,gnval,ptemp(i_num))
        
    sign_new = SIGN(1d0,new_value)
    sign_old = SIGN(1d0,old_value)
    sign_two = SIGN(1d0,sign_new*sign_old)

    if (DABS(new_value) <= DABS(old_valuemax)) then
        index_optmax = i_num
        old_valuemax = new_value        
    end if       
    i_num_star = i_num
    if (sign_two.EQ.-1d0) then
        index_opt = i_num
        old_value = new_value
        GOTO 2525
    ENDIF
    old_value = new_value
ENDDO       
   
2525 IF (index_opt.EQ.0) THEN
    bound1=ptemp(MAX(index_optmax-1,1))
    bound2=ptemp(MIN(index_optmax+1,num))
ELSE
    bound1=ptemp(MAX(index_opt-1,1))
    bound2=ptemp(MIN(index_opt,num))        
ENDIF

p0=BrentRoots(ftype,a0,gnval,gnval,bound1,bound2,1d-10,1000,fval,niter2,error)  
hval = (wbar/alpha/p0)**(1d0/(alpha-1d0))


IF (error.EQ.1) THEN 
    p0 = 0d0
    IF (i_num_star.LT.num)  go to 1313
ELSEIF (error.EQ.2) THEN
    PRINT*,'Problem: Counter 2 UNEMP: More Iterations!'
    PAUSE
    STOP
ENDIF

cnval = hval**alpha-gnval
ctval = (p0*omegac/(1d0-omegac))**(1d0/(1d0+mu))*cnval

1013 aux0  = ctval - a0
wval  = wbar  

IF (try.EQ.1) THEN
    tauval=taxrate*(a0+in*p0*hval**alpha)       !!by construction, no need to check tauval.GE.taxbar w/h<1
ELSE
    tauval = p0*gnval-aux0
ENDIF

taxbar = taxrate*(a0+in*p0*hval**alpha)
qbval  = payoff*b0-aux0

IF (DABS(qbmin-qbval).LE.1d-6)THEN
    hres=-33d0
    PRINT*,'Problem w/counter2_u_unemp'
    PAUSE
    RETURN
    PAUSE
ENDIF

DO ibp=1,nb2         !flip upside down (to find highest MINLOC), i.e. highest zero
    q02(ibp) = q1(iax,nb2-ibp+1)
    b2(ibp)  = bgrid2(nb2-ibp+1)
ENDDO

aux2   = MINLOC(q02*(b2-(1d0-lambda)*b0))        
iqbmin = aux2(1)
iqbmin = nb2-iqbmin+1

IF ((qbmin-qbval.GT.1d-6)) THEN 
    bout = 1
ELSE
    iaux   = nb2

    DO WHILE ((iaux.GT.iqbmin).AND.(q1(iax,iaux)*(bgrid2(iaux)-(1d0-lambda)*b0).GT.qbval))        !find highest index with min value
        iaux = iaux-1
    ENDDO

    bout = 0
    bpval  = LINEAR_INT(q1(iax,iaux:iaux+1)*(bgrid2(iaux:iaux+1)-(1d0-lambda)*b_initial2),bgrid2(iaux:iaux+1),qbval) 
    qval   = qbval/(bpval-(1d0-lambda)*b0)
    
    IF (qval.LT.0d0) THEN
        PRINT*,'Problem price q = 0 in counter2 FE'
        PAUSE
        STOP
    ENDIF
ENDIF

flag00 = ((p0.LE.0d0).OR.(hval.LT.0d0).OR.(hval.GT.1d0).OR.(ctval.LE.0d0).OR.(cnval.LE.0d0).OR.(bout.EQ.1).OR.(tauval.GT.taxbar))    

IF (.NOT. flag00) THEN
    feas = 1d0
ELSE
    feas = 0d0
ENDIF

!negative penalty
penalty = 1d3*((p0.LE.0d0)*(0d0-p0))
penalty = penalty+1d9*((bpval.LT.bmin)*(bmin-bpval))+1d2*(bpval.LT.bmin)
penalty = penalty +1d3*((hval.GT.1d0)*(hval-1d0)+(hval.LT.0d0)*(0d0-hval))

IF (mu .NE. 0d0) THEN
    cval = (omegac*ctval**(-mu)+ (1d0-omegac)*cnval**(-mu))**(-1d0/mu)
ELSE
    cval = (ctval**(omegac))*(cnval**(1d0-omegac))
ENDIF

IF ((cnval.LE.1d-6).OR.(ctval.LE.1d-6)) cval = 0d0

aahat  = (hval+(1d0-hval)*alphac)**(-1d0)
aa1    = hval+(1d0-hval)*(alphac**(1d0-sigma))

IF (sigma.EQ.1d0) THEN
    utilc = (aahat**(1d0-sigma))*aa1*DLOG(MAX(cval,1d-8))
ELSE    
    utilc = (aahat**(1d0-sigma))*aa1*MAX(cval,1d-8)**(1d0-sigma)/(1d0-sigma)
ENDIF
  
IF (sigmag.EQ.1d0) THEN
    utilg = DLOG(MAX(gnval,1d-8))
ELSE
    utilg = MAX(gnval,1d-8)**(1d0-sigmag)/(1d0-sigmag)
ENDIF
    
uval = (1d0-chi)*utilc+chi*utilg
uval = feas*uval + (1d0-feas)*huge_n + penalty

ucost = defaultgrid(i_default_global2)*udef(a_initial2)
uval = uval + ucost

exp_v_next = beta*cont_value(bpval,a_initial2,i_default_global2)
uval = uval + exp_v_next

IF (try.EQ.1) THEN
    !h<1
    IF (root.EQ.0) THEN
        IF (.NOT.flag00) THEN
            hvalx   = hval
            uvalx   = uval
            bpvalx  = bpval
            ctvalx  = ctval
         ELSE
            uvalx = -1d10
        ENDIF
        root=1
    ELSEIF (root.GT.0) THEN
        IF (uval.GT.uvalx) THEN
            hvalx   = hval
            uvalx   = uval
            bpvalx  = bpval
            ctvalx  = ctval
        ENDIF
        root=root+1
    ENDIF
    
    IF (i_num_star.LT.num)     GO TO 1313
ENDIF


IF (try.EQ.1) THEN
    try=2
    
    !h=1
    hval  = 1d0
    p0    = wbar/alpha
    gnval  = gn_initial2
    
    cnval = 1d0-gnval
    ctval = (p0*omegac/(1d0-omegac))**(1d0/(1d0+mu))*cnval
    
    GOTO 1013
ELSEIF (try.EQ.2) THEN   
    IF (uvalx.GT.uval) THEN
        hval   = hvalx
        uval   = uvalx
        bpval  = bpvalx
        ctval  = ctvalx
    ENDIF
ENDIF

bpres= bpval
hres = hval
gnres=gnval
wres = uval

!x_global2  = hres
!gn_global2 = gnres
!counter2_u_unemp = wres
!b_next_global2 = bpres

end SUBROUTINE counter2_u_unemp
!****************************************************************
REAL(8) FUNCTION counter2_u_fe_unc(gn_initial2,iax,a_initial2,b_initial2,i_default_global2,x)

USE GLOBAL

IMPLICIT NONE

REAL(8),INTENT(IN)::a_initial2,b_initial2,gn_initial2,x
INTEGER,INTENT(IN)::i_default_global2,iax
REAL(8)::a0,b0,wres,utilc,utilg,uval
REAL(8)::bpval,wval,tauval,aux0
REAL(8)::p0,h0,q,b
REAL(8)::ctval,cnval,gnval,hval,cval,penalty,fval,ctmax
INTEGER(4)::i,j,jj,feas ,error,niter2,try,numval,ftype
LOGICAL::flag00
REAL(8)::x_glob,taxbar
REAL(8)::qres,ucost,udef,exp_v_next,cont_value,q_fun,old_value,new_value,ftol
INTEGER(4)::im,ibp,iaux,iqbmin,bout,ibx,i_num,index_opt
EXTERNAL cont_value,q_fun
INTEGER,PARAMETER::num=300
REAL(8)::htemp(num)
REAL(8)::den

b0 = b_initial2
a0 = DEXP(a_initial2)
gnval  = gn_initial2
b = x
q = q_fun(b,a_initial2)

aux0 = -q*(b-(1d0-lambda)*b0) + payoff*b0            !p*gn - T
aux0 = aux0*(1d0-defaultgrid(i_default_global2))

ctval  = a0+aux0
hval   = 1d0
cnval  = 1d0-gnval

p0     = (1d0-omegac)/omegac*(ctval/cnval)**(1d0+mu)
tauval = p0*gnval-aux0
taxbar = taxrate*(a0+in*p0)
wval   = alpha*p0

flag00 = ((p0.LE.0d0).OR.(ctval.LE.0d0).OR.(wval.LT.wbar))    

IF (.NOT. flag00) THEN
    feas = 1d0
ELSE
    feas = 0d0
ENDIF

!negative penalty
penalty = 1d3*((p0.LE.0d0)*(0d0-p0)+(tauval.GT.taxbar)*(tauval-taxbar)+(wbar.GT.wval)*(wbar-wval))
penalty = penalty+1d9*((b.LT.bmin)*(bmin-b))+1d2*(b.LT.bmin)

IF (mu .NE. 0d0) THEN
    cval = (omegac*ctval**(-mu)+ (1d0-omegac)*cnval**(-mu))**(-1d0/mu)
ELSE
    cval = (ctval**(omegac))*(cnval**(1d0-omegac))
ENDIF

IF ((cnval.LE.1d-6).OR.(ctval.LE.1d-6)) cval = 0d0

IF (sigma.EQ.1d0) THEN
    utilc = DLOG(MAX(cval,1d-8))
ELSE    
    utilc = MAX(cval,1d-8)**(1d0-sigma)/(1d0-sigma)
ENDIF
    
IF (sigmag.EQ.1d0) THEN
    utilg = DLOG(MAX(gnval,1d-8))
ELSE
    utilg = MAX(gnval,1d-8)**(1d0-sigmag)/(1d0-sigmag)
ENDIF
    
uval = (1d0-chi)*utilc+chi*utilg
uval = feas*uval + (1d0-feas)*huge_n + penalty

ucost = defaultgrid(i_default_global2)*udef(a_initial2)
uval = uval + ucost 

exp_v_next = beta*cont_value(b,a_initial2,i_default_global2)

uval = uval + exp_v_next 

IF (choice_glob.EQ.1) x_global=1d0

wres = uval
counter2_u_fe_unc = wres
RETURN

END FUNCTION counter2_u_fe_unc
!****************************************************************
REAL(8) FUNCTION counter2_u_unemp_unc(gn_initial2,iax,a_initial2,b_initial2,i_default_global2,x)

USE GLOBAL

IMPLICIT NONE

REAL(8),INTENT(IN)::a_initial2,b_initial2,gn_initial2,x
INTEGER,INTENT(IN)::i_default_global2,iax
INTEGER,PARAMETER::num=300
REAL(8)::pres,cnres,ctres,cres,a0,b0,wres
REAL(8)::taures,ares,wwres,utilc,utilg,uval
REAL(8)::pp1,bpval,wval,tauval,aux0
REAL(8)::p0x,bpvalx,wvalx,tauvalx,uvalx,hvalx,hf1,ctvalx,cnvalx
REAL(8)::p0,h0,q,b,hlow,hhgh,hzero,aa1,aahat
REAL(8)::ctval,cnval,gnval,hval,cval,penalty,fval,ctmax
INTEGER(4)::i,j,jj,feas ,error,niter2,try,numval,ftype
LOGICAL::flag00
REAL(8)::q02(nb2),b2(nb2),qval,qbval,qbmin,qbmax,aux2(1),fvalue
REAL(8)::bound1,bound2,x_glob,taxbar
REAL(8)::qvalx,qres,ucost,udef,exp_v_next,cont_value,q_fun,htemp(num),old_value,new_value,ftol
INTEGER(4)::im,ibp,iaux,iqbmin,bout,ibx,i_num
INTEGER(4)::index_opt,i_num_star,i_num_star2,index_optmax,root
real(8)::sign_new,sign_old,sign_two,ctemp2,difc,old_valuemax
EXTERNAL cont_value,q_fun

i_num_star2 = 0
root=0

b0 = b_initial2
a0 = DEXP(a_initial2)
gnval  = gn_initial2

b = x
q = q_fun(b,a_initial2)

aux0 = -q*(b-(1d0-lambda)*b0) + payoff*b0            !p*gn - T
aux0 = aux0*(1d0-defaultgrid(i_default_global2))

ctval  = a0+aux0

try = 1

ftype=36
hlow = (gnval+1d-4)**(1d0/alpha)
hhgh = 1d0

i_num_star=0
index_opt = 0
index_optmax = 0
difc = (hhgh-hlow)/(REAL(num)-1d0)

old_value = AFunction(ftype,a0,gnval,ctval,hlow)
old_valuemax = 10d+33 
       
1414 DO i_num=1+i_num_star,num
    ftype=36
    htemp(i_num) = hlow+(hhgh-hlow)*(REAL(i_num)-1d0)/(REAL(num)-1d0)     
    new_value = AFunction(ftype,a0,gnval,ctval,htemp(i_num))
        
    sign_new = SIGN(1d0,new_value)
    sign_old = SIGN(1d0,old_value)
    sign_two = SIGN(1d0,sign_new*sign_old)

    if (DABS(new_value) <= DABS(old_valuemax)) then
        index_optmax = i_num
        old_valuemax = new_value        
    end if      
        
    i_num_star = i_num
     if (sign_two.EQ.-1d0) then
        index_opt = i_num
        old_value = new_value
        GOTO 3535
    ENDIF
    old_value = new_value
ENDDO       
 

3535 IF (index_opt.EQ.0) THEN
    bound1=htemp(MAX(index_optmax-1,1))
    bound2=htemp(MIN(index_optmax+1,num))
ELSE
    bound1=htemp(MAX(index_opt-1,1))
    bound2=htemp(MIN(index_opt,num))        
ENDIF

hzero=BrentRoots(ftype,a0,gnval,ctval,bound1,bound2,1d-8,1000,fval,niter2,error)  
hval = hzero

IF (error.EQ.1) THEN 
    hval = 0d0
    IF (i_num_star.LT.num)  go to 1414
ELSEIF (error.EQ.2) THEN
    PRINT*,'Problem: Counter 2 UNEMP: More Iterations!'
    PAUSE
    STOP
ENDIF

cnval = hval**alpha-gnval
p0    = (1d0-omegac)/omegac*(ctval/cnval)**(1d0+mu)

1113 wval   = wbar  
tauval = p0*gnval-aux0
taxbar = taxrate*(a0+in*p0*hval**alpha)

flag00 = ((p0.LE.0d0).OR.(hval.LT.0d0).OR.(hval.GT.1d0).OR.(ctval.LE.0d0).OR.(cnval.LE.0d0).OR.(bout.EQ.1).OR.(tauval.GT.taxbar))    

IF (.NOT. flag00) THEN
    feas = 1d0
ELSE
    feas = 0d0
ENDIF

!negative penalty
penalty = 1d3*((p0.LE.0d0)*(0d0-p0))
penalty = penalty+1d9*((b.LT.bmin)*(bmin-b))+1d2*(b.LT.bmin)
penalty = penalty +1d3*((hval.GT.1d0)*(hval-1d0)+(hval.LT.0d0)*(0d0-hval))

IF (mu .NE. 0d0) THEN
    cval = (omegac*ctval**(-mu)+ (1d0-omegac)*cnval**(-mu))**(-1d0/mu)
ELSE
    cval = (ctval**(omegac))*(cnval**(1d0-omegac))
ENDIF

aahat  = (hval+(1d0-hval)*alphac)**(-1d0)
aa1    = hval+(1d0-hval)*(alphac**(1d0-sigma))

IF (sigma.EQ.1d0) THEN
    utilc = (aahat**(1d0-sigma))*aa1*DLOG(MAX(cval,1d-8))
ELSE    
    utilc = (aahat**(1d0-sigma))*aa1*MAX(cval,1d-8)**(1d0-sigma)/(1d0-sigma)
ENDIF    

IF (sigmag.EQ.1d0) THEN
    utilg = DLOG(MAX(gnval,1d-8))
ELSE
    utilg = MAX(gnval,1d-8)**(1d0-sigmag)/(1d0-sigmag)
ENDIF
    
uval = (1d0-chi)*utilc+chi*utilg
uval = feas*uval + (1d0-feas)*huge_n + penalty

ucost = defaultgrid(i_default_global2)*udef(a_initial2)
uval = uval + ucost 

exp_v_next = beta*cont_value(b,a_initial2,i_default_global2)

uval = uval + exp_v_next 

IF (try.EQ.1) THEN
    !h<1
    IF (root.EQ.0) THEN
        IF (.NOT.flag00) THEN
            hvalx   = hval
            uvalx   = uval
            bpvalx  = bpval
            ctvalx  = ctval
         ELSE
            uvalx = -1d10
        ENDIF
        root=1
    ELSEIF (root.GT.0) THEN
        IF (.NOT.flag00) THEN
            IF (uval.GT.uvalx) THEN
                hvalx   = hval
                uvalx   = uval
                bpvalx  = bpval
                ctvalx  = ctval
            ENDIF
        ENDIF
        
        root=root+1
    ENDIF
    
    IF (i_num_star.LT.num)     GO TO 1414
ENDIF


IF (try.EQ.1) THEN
    try=2
    
    !h=1
    hval  = 1d0
    cnval = 1d0-gnval
    p0    = (1d0-omegac)/omegac*(ctval/cnval)**(1d0+mu)
    
    bout  = 1
    IF (DABS(p0-wbar/alpha).LT.1d-10) bout=0
    
    GOTO 1113
ELSEIF (try.EQ.2) THEN
    IF (uvalx.GT.uval) THEN
        hval   = hvalx
        uval   = uvalx
        bpval  = bpvalx
        ctval  = ctvalx
    ENDIF
ENDIF

IF (choice_glob.EQ.1) x_global=hval

wres = uval
counter2_u_unemp_unc = wres
RETURN

!x_global2  = hres
!gn_global2 = gnres
!counter2_u_unemp = wres
!b_next_global2 = bpres

END FUNCTION counter2_u_unemp_unc
!****************************************************************
DOUBLE PRECISION function cont_value(b,a_initial,i_default_global)

USE GLOBAL

IMPLICIT NONE

DOUBLE PRECISION,INTENT(IN)::b,a_initial
INTEGER,INTENT(IN)::i_default_global
INTEGER :: other_type, i,j,h, t, d, NINTV
REAL(8)::u_fun,q_fun,acum,exp_v_next,value_next,a_next,a_fun,scalar
REAL(8)::a_threshold,a_max1,a_min1,a_next1,acum2,acum1,g,q,DNORDF,v0_fun,vaut_fun
REAL(8)::acum_int,acum_int_lo,acum_int_hi,borrowing, interpolate,cdf,cdf1,cdf_threshold, DNORIN,objective_u_unemp
REAL(8)::exp_v_next_aut,acum3,objective_u_fe
EXTERNAL u_fun,q_fun,DNORDF,a_fun,interpolate,v0_fun,vaut_fun,DNORIN,udef,objective_u_unemp,objective_u_fe

REAL(8)::value1,value2 ,u_value,ucost,udef

IF (psi1.LT.500d0) THEN
    a_threshold   = a_fun(b,a_initial)
ELSE
    a_threshold   = -1d4
ENDIF

cdf_threshold = DNORDF((a_threshold - rhoa*a_initial - (1-rhoa)* mua)/sigmaa)

exp_v_next = 0d0
exp_v_next_aut = 0d0
scalar = 1d0 / SQRT(pi_number)

a_min1 = (1-rhoa)*mua + rhoa* a_initial - 4*sigmaa
a_max1 = (1-rhoa)*mua + rhoa* a_initial + 4*sigmaa

NINTV = ncdf - 1
acum1=0d+0
acum2=0d+0

   if (a_threshold < a_min1) then

       do i=1,nquad_v
          cdf = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdfmin) + cdfmin
          a_next  = (1-rhoa)*mua + rhoa* a_initial + DNORIN(cdf) * sigmaa
          value_next = v0_fun(b,a_next)
          acum1 = acum1 + 0.5*(cdfmax-cdfmin)*quad_w_v(i) * value_next
       end do

   elseif(a_threshold > a_max1) then

       do i=1,neps
          cdf = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdfmin) + cdfmin
          a_next  = (1-rhoa)*mua + rhoa* a_initial + DNORIN(cdf) * sigmaa
          value_next = vaut_fun(a_next)
          acum1 = acum1 + 0.5*(cdfmax - cdfmin)*quad_w_v(i) * value_next
       end do

   else
      do i=1,neps

        cdf = 0.5*(quad_x_v(i) + 1)*(cdf_threshold - cdfmin) + cdfmin
        a_next  = (1-rhoa)*mua + rhoa* a_initial + DNORIN(cdf) * sigmaa
        value_next = vaut_fun(a_next)
        acum1 = acum1 + 0.5*(cdf_threshold - cdfmin)*quad_w_v(i) * value_next

!        compute integral for a > a_threshold
        cdf1 = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdf_threshold) + cdf_threshold
        a_next1 = (1-rhoa)*mua + rhoa* a_initial + DNORIN(cdf1) * sigmaa 
        value_next = v0_fun(b,a_next1)
        acum2 = acum2 + 0.5*(cdfmax - cdf_threshold)*quad_w_v(i) * value_next
      end do
   end if

exp_v_next = (acum1+acum2)/(cdfmax - cdfmin)

IF (i_default_global.EQ.2) THEN
    acum3=0d+0
        do i=1,neps
            cdf = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdfmin) + cdfmin
            a_next  = (1-rhoa)*mua + rhoa* a_initial + DNORIN(cdf) * sigmaa
            value_next = vaut_fun(a_next)
            acum3 = acum3 + 0.5*(cdfmax - cdfmin)*quad_w_v(i) * value_next
        end do
    
        exp_v_next_aut = acum3/(cdfmax - cdfmin)
        exp_v_next = phi*exp_v_next + (1d0-phi)*exp_v_next_aut
ENDIF

cont_value = exp_v_next

END FUNCTION
!********************************************************************************************
subroutine optimize_fe_unc(b_next,v_value,gn_initial2,iax,a_initial2,b_initial2,i_default_global2)
USE GLOBAL

implicit none

REAL(8),INTENT(IN)::a_initial2,b_initial2,gn_initial2
INTEGER,INTENT(IN)::i_default_global2,iax
integer :: MAXFN, i,t, d, index_opt, index_vector(1), num, num_tirar
parameter (num=750, num_tirar = 500)
DOUBLE PRECISION :: b_next, v_value, counter2_u_fe_unc, new_value, q_fun, old_value, b, u_fun,q,vector(nb),&
                    b_next_grid(num), STEP, BOUND, XACC, b_next_guess
real(8)::b_next_inf,b_next_sup 
EXTERNAL counter2_u_fe_unc, q_fun, u_fun, DUVMIF

b_next_inf = bmin
b_next_sup = bmax

old_value = 10d+33
do i=1,num
   b_next_grid(i) = b_next_inf + (b_next_sup - b_next_inf) * (i-1)/ (num - 1)
   new_value = objective_fun4(b_next_grid(i))
   if (new_value <= old_value) then
      index_opt = i
      b_next_guess = b_next_grid(i)
      old_value = new_value
  end if
end do

STEP = (b_next_grid(2) - b_next_grid(1))*0.5
BOUND = 10*EXP(amax)
XACC = 1d-6
MAXFN = 1000
   
IF (old_value.GE.1d5) THEN
    b_next = 0d0
    v_value = -old_value
ELSE
    call DUVMIF(objective_fun4, b_next_guess, STEP, BOUND, XACC, MAXFN, b_next)
    v_value = -objective_fun4(b_next)     
end if

CONTAINS

REAL(8) FUNCTION objective_fun4(x)

REAL(8),INTENT(IN)::x

objective_fun4 = -counter2_u_fe_unc(gn_initial2,iax,a_initial2,b_initial2,i_default_global2,x)
END FUNCTION

end subroutine
!*****************************************************************
subroutine optimize_unc(b_next,v_value,gn_initial2,iax,a_initial2,b_initial2,i_default_global2)
USE GLOBAL

implicit none
 
REAL(8),INTENT(IN)::a_initial2,b_initial2,gn_initial2
INTEGER,INTENT(IN)::i_default_global2,iax
integer :: MAXFN, i,t, d, index_opt, index_vector(1), num, num_tirar
parameter (num=750, num_tirar = 500)
DOUBLE PRECISION :: b_next, v_value, counter2_u_unemp_unc, new_value, q_fun, old_value, b, u_fun,q,vector(nb),&
                    b_next_grid(num), STEP, BOUND, XACC, b_next_guess
real(8)::b_next_inf,b_next_sup 
EXTERNAL counter2_u_unemp_unc, q_fun, u_fun, DUVMIF

b_next_inf = bmin
b_next_sup = bmax

old_value = 10d+33

do i=1,num
    b_next_grid(i) = b_next_inf + (b_next_sup - b_next_inf) * (i-1)/ (num - 1)
    new_value = objective_fun3(b_next_grid(i))
    if (new_value <= old_value) then
      index_opt = i
      b_next_guess = b_next_grid(i)
      old_value = new_value
    end if
end do

STEP = (b_next_grid(2) - b_next_grid(1))*0.5
BOUND = 10*EXP(amax)
XACC = 1d-6
MAXFN = 1000
   
IF (old_value.GE.1d5) THEN
    b_next = 0d0
    v_value = -old_value
ELSE
    call DUVMIF(objective_fun3, b_next_guess, STEP, BOUND, XACC, MAXFN, b_next)
    v_value = -objective_fun3(b_next)    
end if

CONTAINS

REAL(8) FUNCTION objective_fun3(x)

REAL(8),INTENT(IN)::x

objective_fun3 = -counter2_u_unemp_unc(gn_initial2,iax,a_initial2,b_initial2,i_default_global2,x)
END FUNCTION

end subroutine
!******************************************************************************************
SUBROUTINE counter1(iax,ibx)
! Counterfactual exercise fixing (A,b) 
! Compute prices, policy and util fcns for alternative values of gN financed through debt and taxes
! istate = 0 repayment today, =1 default today

USE GLOBAL

IMPLICIT NONE

INTEGER(4)::iax,ibx
REAL(8)::gnn2(nn),aux00,aux11,aux12,gnoptimal,p0,b_next,gn_value,h_value,v_value
REAL(8)::qoff(nn),hoff(nn),gnoff(nn),ctoff(nn),woff(nn),bpoff(nn),tauoff(nn),wwoff(nn)
REAL(8)::q_fun,q,value1,value2,value3,value4,u_value 
REAL(8)::b_next1,b_next2,b_next3,b_next4,h_value1,h_value3,h_value4,gn_value1,h_value2,gn_value2,p00
REAL(8)::a_initial,b_initial,gn_initial,a0, aux0
INTEGER::i_default_global 

EXTERNAL::q_fun,counter2_u_fe,counter2_u_unemp
INTEGER::nbaux,ign,i_b0(1)

nbaux = nb2

!construct grid for gn
aux11 = MIN(gn1(iax,ibx),gnaut1(iax))
aux12 = 1d0-MAX(gn1(iax,ibx),gnaut1(iax))
    
aux00  = 0.65d0*MIN(aux11,aux12)
aux00  = 2d0*aux00/DBLE(nn-1)
    
gnoptimal = gn1(iax,ibx)   
gnn2(1)   = gnoptimal -aux00*(nn-1)/2d0
 
DO ign=2,nn
    gnn2(ign)=  gnn2(ign-1)+aux00   
ENDDO

!set state vector
a_initial = agrid2(iax)
b_initial = bgrid2(ibx)

i_default_global = 1 !repays

a0 = DEXP(a_initial)

choice_glob=1

DO ign=1,nn
    gn_initial = gnn2(ign)

    CALL counter2_u_fe(gn_initial,iax,a_initial,b_initial,i_default_global,b_next1,h_value1,gn_value1,value1)
    CALL counter2_u_unemp(gn_initial,iax,a_initial,b_initial,i_default_global,b_next2,h_value2,gn_value2,value2)
    call optimize_unc(b_next3,value3,gn_initial,iax,a_initial,b_initial,i_default_global)
    h_value3= x_global
    call optimize_fe_unc(b_next4,value4,gn_initial,iax,a_initial,b_initial,i_default_global)
    v_value = MAX(value1,MAX(value4,MAX(value2,value3)))
   
    
    IF (value1.EQ.v_value) THEN
        !reg_global = 1
        h_value    = h_value1  
        gn_value   = gn_initial
        b_next     = b_next1    
        
        q = q_fun(b_next,a_initial)
        aux0 = -q*(b_next-(1d0-lambda)*b_initial) + payoff*b_initial           !p*gn - T
                
        p0 = (1d0-omegac)/omegac*((a0+aux0)/(1d0-gn_value))**(1d0+mu)
       
    ELSEIF (value2.EQ.v_value) THEN
        
        !reg_global = 2
        h_value    = h_value2
        gn_value   = gn_initial
        b_next     = b_next2
        p0         = wbar/alpha*h_value**(1d0-alpha)
            
        q = q_fun(b_next,a_initial)        
        aux0 = -q*(b_next-(1d0-lambda)*b_initial) + payoff*b_initial           !p*gn - T
    
    ELSEIF (value3.EQ.v_value) THEN

        h_value    = h_value3
        gn_value   = gn_initial
        b_next     = b_next3
                
        q = q_fun(b_next,a_initial)        
        aux0 = -q*(b_next-(1d0-lambda)*b_initial) + payoff*b_initial           !p*gn - T
        p0 = (1d0-omegac)/omegac*((a0+aux0)/(h_value**alpha-gn_value))**(1d0+mu)

    ELSEIF (value4.EQ.v_value) THEN
        !reg_global = 1
        h_value    = 1d0  
        gn_value   = gn_initial
        b_next     = b_next4    
        
        q = q_fun(b_next,a_initial)
        aux0 = -q*(b_next-(1d0-lambda)*b_initial) + payoff*b_initial           !p*gn - T
                
        p0 = (1d0-omegac)/omegac*((a0+aux0)/(1d0-gn_value))**(1d0+mu)
       
    ENDIF 
    
    qoff(ign)  = q
    hoff(ign)  = h_value
    IF (v_value.LE.-1d3) hoff(ign)=0d0

    gnoff(ign) = gn_value
    ctoff(ign) = a0+aux0
    woff(ign)  = v_value
    
    bpoff(ign) = b_next
    
    p0 = (1d0-omegac)/omegac*(ctoff(ign)/(h_value**alpha-gn_value))**(1d0+mu)
    tauoff(ign)= p0*gn_value-aux0        
    wwoff(ign) = MIN(10d0,alpha*p0*h_value**(alpha-1d0))
 
    IF (h_value.EQ.0d0) wwoff(ign)=0d0
ENDDO


OPEN(1,FILE='output//counter//rep//hoff.txt')
OPEN(2,FILE='output//counter//rep//gnoff.txt')
OPEN(3,FILE='output//counter//rep//ctoff.txt')
OPEN(4,FILE='output//counter//rep//woff.txt')
OPEN(66,FILE='output//counter//rep//wwoff.txt')
OPEN(7,FILE='output//counter//rep//bpoff.txt')
OPEN(8,FILE='output//counter//rep//tauoff.txt')
OPEN(9,FILE='output//counter//rep//g.txt')
OPEN(10,FILE='output//counter//rep//qoff.txt')
OPEN(11,FILE='output//counter//rep//waut.txt')      !relevant to determine default

DO ign=1,nn
    WRITE(1,*) hoff(ign)
    WRITE(2,*) gnoff(ign)
    WRITE(3,*) ctoff(ign)
    WRITE(4,*) woff(ign)
    WRITE(66,*) wwoff(ign)
    WRITE(7,*) bpoff(ign)
    WRITE(8,*) tauoff(ign)
    WRITE(9,*) gnoff(ign)    
    WRITE(10,*) qoff(ign)
ENDDO

WRITE(11,*) vaut1(iax)

CLOSE(1)
CLOSE(2)
CLOSE(3)
CLOSE(4)
CLOSE(66)
CLOSE(7)
CLOSE(8)
CLOSE(9)
CLOSE(10)
CLOSE(11)


END SUBROUTINE
!***************************************************************************
SUBROUTINE fiscal_multi

USE GLOBAL

IMPLICIT NONE

INTEGER(4)::i_a,i_b
REAL(8)::gnn2(3),aux00,aux11,aux12,gnoptimal,p0,b_next,gn_value,h_value,v_value
REAL(8)::q_fun,q,value1,value2,value3,value4,u_value 
REAL(8)::b_next1,b_next3,h_value1,h_value3,gn_value1,b_next2,b_next4,h_value2,gn_value2,p00
REAL(8)::a_initial,b_initial,gn_initial,a0, aux0
INTEGER::i_default_global 
REAL(8)::ctval,cnval

EXTERNAL::q_fun,counter2_u_fe,counter2_u_unemp
INTEGER::nbaux,ign,i_b0(1)

REAL(8)::vargdp,varvgn,vargn,varyn,varw,varr,varq,varbp  
REAL(8)::fm(na2,nb2),valf(na2,nb2),valfaut(na2),rchange(na2,nb2),derivb(na2,nb2)
REAL(8)::woff(3),hoff(3),poff(3),gnoff(3),gdpoff(3),roff(3),qoff(3),bpoff(3)
REAL(8)::new_value,max_case,derb
REAL(8)::qoptimal(na2,nb2),qgN2(na2,nb2),bpoptimal(na2,nb2),bpgN2(na2,nb2)     
    
INTEGER(4)::index_bfeas22,index_opt,index_a,index_b

REAL(8)::fm1,fmn1,valf1,gn_basex(na,501) ,bgridx(501)
index_bfeas4 = 0d0
        
!Fiscal mutipliers computed using pNss
!debt indices:

!default 
!bss=107
!meanb=122

!flex wages:
!bss=156

index_a=11
index_b=107
i_b0=MINLOC((bgrid2(bp1(index_a,:))-bgrid2)**2d0)
index_b=i_b0(1)


!to compute spline for feasible debt gridpoints
do i_a = 1,na
    new_value=-1d5
    i_b=0
    index_opt = 0
        
    do while (new_value <= -1d3)
        index_opt = i_b
        IF (i_b.EQ.nb) EXIT
        i_b=i_b+1
        new_value = v0_matrix(i_b,i_a)
    end do
    index_bfeas4(i_a)=index_opt        
ENDDO

index_bfeas22=MAXVAL(index_bfeas4)

do i_b = 1,nb2  
    do i_a = 1,na2
    IF ((vcon1(i_a,i_b).LE.-1d2).OR.(def1(i_a,i_b).EQ.0))  THEN       
        fm(i_a,i_b)  = -1d0
        valf(i_a,i_b)= 0d0
    ELSE        
         
        a_initial = agrid2(i_a)
        b_initial = bgrid2(i_b) 

        i_default_global = 1 !repays                 

        aux11 = gn1(i_a,i_b)  
        aux12 = 1d0-aux11
                
        gnoptimal = aux11   
        gnn2(1)   = MAX(1d-5,aux11-0.001d0)
        gnn2(2)   = aux11+0.001d0        
        gnn2(3)   = aux11        

        a0 = DEXP(a_initial)

        choice_glob=1
        DO ign=2,3      
            gn_initial = gnn2(ign)

            CALL counter2_u_fe(gn_initial,i_a,a_initial,b_initial,i_default_global,b_next1,h_value1,gn_value1,value1)
            CALL counter2_u_unemp(gn_initial,i_a,a_initial,b_initial,i_default_global,b_next2,h_value2,gn_value2,value2)            
            call optimize_fe_unc(b_next4,value4,gn_initial,i_a,a_initial,b_initial,i_default_global)
            call optimize_unc(b_next3,value3,gn_initial,i_a,a_initial,b_initial,i_default_global)    
            v_value = MAX(value1,MAX(value4,MAX(value2,value3)))
            h_value3= x_global

        IF (value1.EQ.v_value) THEN
                h_value    = h_value1  
                gn_value   = gn_initial
                b_next     = b_next1    
        
                q = q_fun(b_next,a_initial)
                aux0 = -q*(b_next-(1d0-lambda)*b_initial) + payoff*b_initial           !p*gn - T
                
                p0 = (1d0-omegac)/omegac*((a0+aux0)/(1d0-gn_value))**(1d0+mu)

            ELSEIF (value4.EQ.v_value) THEN
                h_value    = 1d0  
                gn_value   = gn_initial
                b_next     = b_next4    
        
                q = q_fun(b_next,a_initial)
                aux0 = -q*(b_next-(1d0-lambda)*b_initial) + payoff*b_initial           !p*gn - T
                
                p0 = (1d0-omegac)/omegac*((a0+aux0)/(1d0-gn_value))**(1d0+mu)
       
            ELSEIF (value2.EQ.v_value) THEN

                h_value    = h_value2
                gn_value   = gn_initial
                b_next     = b_next2
                p0         = wbar/alpha*h_value**(1d0-alpha)
            
                q = q_fun(b_next,a_initial)        
                aux0 = -q*(b_next-(1d0-lambda)*b_initial) + payoff*b_initial           !p*gn - T
             
            ELSEIF (value3.EQ.v_value) THEN

                h_value    = h_value3
                gn_value   = gn_initial
                b_next     = b_next3
                
                q = q_fun(b_next,a_initial)        
                aux0 = -q*(b_next-(1d0-lambda)*b_initial) + payoff*b_initial           !p*gn - T
                p0 = (1d0-omegac)/omegac*((a0+aux0)/(h_value**alpha-gn_value))**(1d0+mu)

            ENDIF 

            ctval = a0+aux0
            
            hoff(ign)  = h_value
            IF (v_value.LE.-1d3) hoff(ign)=0d0

            gnoff(ign) = gn_value
            woff(ign)  = v_value
    
            poff(ign)    = p1(index_a,index_b)  
            gdpoff(ign)  = a0+poff(ign)*h_value**alpha
            roff(ign)    = aux0 - payoff*b_initial 
            qoff(ign)    = q
            bpoff(ign)   = b_next

        ENDDO

        fm1  = 100d0
        valf1= 100d0

        DO ign=2,2
            IF (gdpoff(ign).GT.0d0) THEN
        
                vargdp = gdpoff(ign)-gdpoff(3)
                varvgn = poff(ign)*gnoff(ign)-poff(3)*gnoff(3)
                vargn  = gnoff(ign)-gnoff(3)
                varyn  = hoff(ign)**alpha-hoff(3)**alpha
                varw   = woff(ign)-woff(3)
                varr   = roff(ign)-roff(3)
                varq   = qoff(ign)-qoff(3)
                varbp  = bpoff(ign)-bpoff(3)
            
                fm1  = MIN(fm1,vargdp/varvgn)
                valf1= MIN(fm1,varw)   
                derb = varq/varbp*(bpoff(3)-(1d0-lambda)*b_initial)
            ENDIF
        ENDDO
    
    fm(i_a,i_b)  = fm1
    valf(i_a,i_b)= valf1
    rchange(i_a,i_b)= varr
    derivb(i_a,i_b) = derb
    qoptimal(i_a,i_b)  = qoff(3)
    qgN2(i_a,i_b)      = qoff(2)
    bpoptimal(i_a,i_b) = bpoff(3)
    bpgN2(i_a,i_b)     = bpoff(2)
    
    ENDIF

    ENDDO

ENDDO
 

OPEN(1,FILE='output//fm_ss.txt')
OPEN(3,FILE='output//fmval_ss.txt')
OPEN(4,FILE='output//rchange_ss.txt')
OPEN(5,FILE='output//derivativeb_ss.txt')
OPEN(6,FILE='output//qoptimal_ss.txt')
OPEN(7,FILE='output//qgn2_ss.txt')
OPEN(8,FILE='output//bpoptimal_ss.txt')
OPEN(9,FILE='output//bpgn2_ss.txt')

do i_b = 1,nb2 
do i_a = 1,na2 
    
    WRITE(1,*) fm(i_a,i_b)  
    WRITE(3,*) valf(i_a,i_b)
    WRITE(4,*) rchange(i_a,i_b)
    WRITE(5,*) derivb(i_a,i_b)
    WRITE(6,*) qoptimal(i_a,i_b) 
    WRITE(7,*) qgN2(i_a,i_b)     
    WRITE(8,*) bpoptimal(i_a,i_b) 
    WRITE(9,*) bpgN2(i_a,i_b)     
ENDDO
ENDDO

CLOSE(1)
CLOSE(3)
CLOSE(4)
CLOSE(5)
CLOSE(6)
CLOSE(7)
CLOSE(8)
CLOSE(9)

PRINT*,'FISCAL MULTIPLIERS UPLOADED'

end SUBROUTINE fiscal_multi

    
    
    
    
    
    
    
    
    
    
    
    
    
